Transport Coefficients of the Yukawa One Component Plasma. 
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Abstract 

We present equilibrium molecular-dynamics computations of the thermal 

conductivity and the two viscosities of the Yukawa one-component plasma. 

The simulations were performed within periodic boundary conditions and 

Ewald sums were implemented for the potentials, the forces, and for all the 

currents which enter the Kubo formulas. For large values of the screening 

parameter, our estimates of the shear viscosity and the thermal conductivity 

are in good agreement with the predictions of the Chapman-Enskog theory. 
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Recently, many numerical studies of the Yukawa one-component plasma (YOCP) - i.e. 
a system made of identical classical point particles of charge q and mass m which are 
embedded in a uniform neutralizing background of volume V and which interact via Yukawa 
pair potentials v{r) = q"^ exp{—ar)/r - have been performed in view of applications for a 
broad variety of systems, including dusty plasmas, inertial confinement fusion dense plasmas, 
jovian planets, brown and white dwarfs, etc. The excess free energy / as well as all the 
thermodynamic properties of the YOCP depend only upon two parameters, namely the 
coupling parameter F = f3q^ /a, where j3 = 1/kT is the inverse temperature and a the ionic 
radius [Airpa^/S = 1, p = N/V number density), and the reduced screening parameter 
a* = aa. In the special case where a* = one recovers the well-known one-component 
plasma (OCP) The other limiting case a* ^ oo is that of a dilute gase for which simple 
approximate schemes can safely be used. The thermodynamic and structural properties 
of the YOCP have been thoroughly studied by means of equilibrium molecular dynamics 
(EMD) simulations within periodic boundary conditions (PBC) and by Monte Carlo 
simulations on the hypersphere Reliable estimates of the free energy f(r,a*) are thus 
available in a wide range of (F, a*) 0-0. 

By contrast, very little is known about the dynamical properties of the model. In view 
of hydrodynamical simulations, precise estimates of the transport coefficients of the YOCP 
are clearly wanted. Attempts to compute the shear viscosity rj by means of non-equilibrium 
molecular dynamics (NEMD) simulations were discussed recently in the literature ||^. In 
this letter we present (EMD) computations not only of f], but also of the bulk viscosity ^ 
and the thermal conductivity A. It turns out that our results for rj differ significantly from 
those of ref . , a puzzling point which will be discussed further on. 

As it is well known, the three transport coefficients r/, ^, and A are given by the Kubo 
formulas : 




(la) 
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^=^E/ {^o.a{t)appm dt , (lb) 



A = [ {Je{t) ■ Jem dt . (IC) 







In Eqs. (1) aai3 denotes the Fourier transform of one of the cartesian components of the 
pressure tensor at = and Je is the A; = component of the Fourier transform of the 
energy current. 

Our simulations were performed in a cube of side L with PBC conditions and we took 
an exphcit account of the periodicity of the system by making use of Ewald sums. We 
have shown elsewhere that the PBC expression of the Yukawa pair potential reads, up to an 
additional constant, as 0: 

47rg2 exp(-(k^ + a^) /46^) , ^ erfc((5||f1| + ea/2(5) exp(ea||r1 

VPBcir) = —f- J2 7^ Lexp{ik-r) + q^ J2 



(2) 

where the sum in the r.h.s runs over the vectors k of the reciprocal lattice. The parameter 
6 is chosen in such a way that the contributions to vpbc{^ in the direct space reduce to a 
single term ( i.e. the second term of the r.h.s. of Eq. (H)) and that the cut-off ko on the 
vectors k is not too large. The optimal choice, which ensures a relative precision of the order 
of ~ 10~^ on vpBC'i'f^) the points r inside the simulation cell is 5 x L ~ 5.6 

The Ewald expressions for the pressure tensor aa/3 and the energy current Je can be 
obtained by generalizing the pioneer work of Bernu and Vieillefosse on the OCP [|]. The 
details of the derivation will be given elsewhere, and we just quote here, as an example, the 
resulting expression for the k = Fourier transform of the the pressure tensor : 

^a,/3 = </3 + ^1/3 + . (3a) 



N 

^a,l3 = ^Y1 , (3b) 

i=l 
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(i erfc(5||r|| + ea;/2(5) exp(ea;r) 



\ dr 2r 



(3c) 



^ ^ S „.B.o — ^^^^^ — ' ^ ~ w^^'-'n • 

(3d) 

In our simulations we choose as unit of length the ionic radius a and as unit of time 
r = \/?)UJp^ with ujp = Aixpq^ /m. The calculations were performed in the microcanonical 
ensemble and the trajectories of each of the particles (and all its images) were computed 
by a time-symmetrical integer algorithm This algorithm is symplectic and ensures 

an exact conservation of the total momentum of the system. The time increment At was 
chosen in such a way to ensure a good conservation of the energy (typically At = lO^^r 
leads to fluctuations of ~ 10"'' on the average energy). In most of our simulations = 500, 
but smaller and larger systems were also considered in order to study finite size effects 
on the transport coefficients. Typically 5 x 10^ time steps were generated after a careful 
equilibration of the system. Each run was divided into statistically independent blocks of 
~ 5 X 10'^ time steps, i.e. much larger than the correlation time. The reported errors on the 
autocorrelation functions and the transport coefficients were obtained by a standard block 
analysis ; they correspond to one standard deviation. As an example of the precision which 
can be obtained for sufficiently long calculations, we display in Fig. |l| the autocorrelation 
of the energy current at (F = 10,0;* = 0.1). The integral of the function over the time 
reaches a well defined plateau which allows for an accurate determination of the thermal 
conductivity. The precision on A and on the other transport coefficients is typically of the 
order of ~ 1% for most of the considered cases. 

Since the thermodynamic states of the YOCP are characterized by two parameters, 
a systematic determination of the transport coefficients in the whole fluid phase requires 
an enormous amount of simulations. We present here only preliminary results for a few 



thermodynamic states; extended results will be given in a forthcoming publication [11 



Our results are summarized in Tables I, II, III. We choose the following units: muppo? 
for the viscosities {rj = mUppa^r]*, ^ = mUppa?^*) , kujppa? for the thermal conductivity 
(A = kuppa^X*). 

In order to check our method we have first considered the case a* = 0.01 and compared 
our results with those of Bernu and Vieillefosse |^ and of Donko et al. for the OOP. The 
former authors have performed EMD simulations of the OCP and give estimates of (r^, ^, A) 
for a few thermodynamic states, while the latter provide extensive NEMD computations of 77 
and A. As far as the shear viscosity is concerned, all the results are in overall good agreement 
except at low F's. However, it must be stressed that, in this regime, Bernu and Vieillefosse 
have considered only relatively small systems of = 128 particles; their results for r] are 
hence probably underestimated due to finite size effects. As seen from table I, the reduced 
bulk viscosity ^* is typically three orders of magnitude smaller than 77* which makes difficult 
its precise determination and entails relatively important statistical errors. Our estimates 
of ^* agree well with those of Bernu et al. at large F's but, as for the shear viscosity, differ 
significantly at low F's. Finally, our results for the thermal conductivity at a* = 0.01 are in 
good agreement with those of Bernu et al. (except for the lowest F's) but are systematically 
higher than those obtained by Donko et al. in their NEMD simulations . 

The recent NEMD simulations of Sanbonmatsu and M. S. Murillo on the YOCP have 
been performed only for large values of the screening parameter ( i.e. for a* = 1, 2, 3, 4). In 
this case, Ewald sums can probably be safely ignored, at least for sufficiently large systems. 
Our estimates of the shear vicosity are compared with those of ref. in table II for a few 
thermodynamic states. The disagreement between the two series of simulations is patent, 
particularly for large values of a* where the results may differ by a factor of ~ 4. In order to 
clarify this point, we have focussed on the case of the large a*'s. In this regime we actually 
deal with a dilute gas of particles interacting via short range potentials. Clearly, in this 
case, the transport coefficients can be computed in the framework of the Chapman-Ensog 
(CE) theory |T^. In the so-called first CE approximation, we have ^* = and 
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where fi(^)(2) is a standard collision integral [|I^]. Note first that = which is compatible 
with the low values of the reported data and the steady decay of ^* with respect to a* for a 
fixed r, as seen from table II. Moreover, it can be shown that the expression of the CE 
shear viscosity of the YOCP can be rewritten as 



«*2 



VCE = ^ X(aT) , (5) 

where X{a*r) is a triple integral that we have computed numerically by Monte Carlo in- 
tegration methods. In Fig. 2 we display the EMD and CE shear viscosities as functions 
of a* for r = 2, 10, 50. The agreement between our EMD results and the predictions of 
the CE theory is obvious for sufficiently large a*'s. More precisely, the CE estimates seem 
to be accurate as soon as the coupling parameter rexp(— a*) ^ 0.35. The CE theory also 
enables the computation of the thermal conductivity \ce = ^C^rjcE/'^ and we found, as in 
the case of the shear viscosity, a perfect agreement between our EMD simulations and the 
CE theoretical predictions in the domain rexp(— a*) ^ 0.35. 

To summarize, our EMD results for the transport coefficients of the YOCP are in good 
agreement with the available data on the OCP in the limit a* — > and also in good 
agreement with the predictions of the CE theory for large values of a*, as it should be, and 
in severe disagreement with the values reported in ref. 

We think that the standard approach used in this work to compute the transport coeffi- 
cients -i.e. EMD simulations plus Ewald sums - is efficient and reliable for the two following 
reasons : 

• the three transport coefficients r], ^, and A can be computed in a single run. By 
contrast, each transport coefficient requires a separate NEMD simulation. 

• thanks to Ewald sums the simulations can be undertaken for any value of a* and they 
require a small number of particles. By contrast NEMD simulations seem to require 
larger system sizes which precludes the use of Ewald sums P,|l2 . 
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We have indeed checked that finite size effects on the transport coeeficients are small as 
long as > 256. For instance, for the state (F = 10, a* = 1), we found for the thermal 
conductivity A* = 0.4138(41), 0.5556(54), 0.5397(69), 0.5372(56) for = 128,256,500,864 
respectively. Therefore systems of ~ 500 are sufficiently large to ensure a reliable estimate 
of the transport coefficients. Some discrepencies between our results in the case a* = 0.01 
and those obtained by Bernu et al. for the OCP with systems made of = 128 particles 
probably originate in finite size effects. 

Finally we have also considered small values of the screening parameter a*, i.e. < a* < 
1. In this case the use of Ewald sums cannot be avoided and some preliminary results are 
diplayed in table III. Calculations are in progrees for other values of (F, a*) and many more 
results will be given together with a fit of all transport coefficients as functions of (F, a*) 

ini- 

We acknowledge D. Gilles, J. Clerouin, and D. Levesque for useful discussions and D. 
Levesque for providing us a MD code of the OCP easily transformed in a MD code for the 
YOCP. 
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FIGURES 

FIG. 1. Solid curve : the autocorrelation function of the energy current {Je{t) ■ ^e(O)) for 
(r = 10, a* = 0.1), dotted curve : cumulative sum. 

FIG. 2. Shear viscosity of the YOCP as a function of a* for various F's. Solid curve : EMD 
results. Dashed curve : CE theory. 
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TABLES 

TABLE L Transport coefficients of the YOCP in the hmit a* 0. The numbers in brackets 
denote the accuracy of the last digits. 



V 



C X 10 



-3 



A 



r Yocp^ 

1 1.16(5) 

2 0.527(7) 



OCP^ OCP= 
1.04(21) 

~ 0.5 



10 0.112(1) 0.085(17) ~ 0.1 
100 0.1874(20) 0.18(3) ~ 0.18 



YOCP^ OCP^ 

4.72(36) 2.6(6) 
3.02(5) 

1.753(24) 1.8(5) 



YOCP^ 
4.24(29) 
1.862(16) 



OCP^ OCP^ 
2.9(6) ~ 2.2 
~ 1.2 



0.5586(70) 0.66(16) 



0.394(7) 0.21(6) 0.843(11) 0.88(17) 



'"EMD results at a* = 0.01 

'^EMD results of Bernu and Vieillefosse Ref. |] for the OCP. 
'^NEMD results of Donko et al. Ref. [HI for the OCP. 



0.40 
0.72 
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TABLE II. Transport coefficients of tfie YOCP for few tfiermodynamic states. For each coeffi- 
cient, first column : present EMD results, second column : Chapmann-Enskog prediction {^CE = 
not reported) , third column (only for r/*) NEMD estimates of ref. [5]. The numbers in brackets 
denote the accuracy of the last digits. 



r = 2 



a 




V 






A 


^ X 10-3 


1 


0.496(12) 


0.439(64) 


0.2340 


2.42(12) 


1.65(24) 


0.834(48) 


2 


0.991(24) 


0.826(48) 


0.2646 


2.89(17) 


3.09(18) 


0.756(14) 


3 


1.282(36) 


1.367(94) 


0.4760 


5.36(29) 


5.13(35) 


0.694(12) 


4 


1.935(36) 


2.055(152) 


0.5496 


7.18(23) 


7.71(57) 


0.447(5) 










r 


= 10 






a 




V 




A 




^ X 10-3 


1 


0.112(3) 


0.047(5) 


0.0526 


0.570(18) 


0.176(18) 


1.282(48) 


2 


0.145(3) 


0.117(10) 


0.0521 


0.644(17) 


0.438(39) 


1.205(48)(14) 


3 


0.198(3) 


0.193(10) 


0.0693 


0.841(18) 


0.726(40) 


1.426(12) 


4 


0.306(4) 


0.288(19) 


0.0870 


1.239(23) 


1.08(7) 


1.255(9) 



TABLE III. Transport coefficients of the YOCP for small values of the screening parameter a. 
The numbers in brackets denote the accuracy of the last digits. 











[: X 10-^ 






A 




a r = 2 


r = 10 


r = 50 


r = 2 


r = 10 


r = 50 


r = 2 


r = 10 


r = 50 


0.2 0.513(5) 


0.1054(12) 


0.1102(7) 


2.78(5) 


1.287(24) 


0.534(7) 


1.716(11) 


0.55(1) 


0.641(1) 


0.4 0.464(5) 


0.1033(12) 


0.1069(7) 


1.439(24) 


1.238(24) 


0.451(5) 


1.96(2) 


0.55(1) 


0.704(1) 


0.6 0.513(5) 


0.1093(17) 


0.1016(6) 


0.967(12) 


0.914(17) 


0.3906(85) 


1.99(2) 


0.518(9) 


0.560(1) 


0.8 0.525(5) 


0.1028(15) 


0.0937(6) 


0.851(5) 


0.788(19) 


0.372(5) 


2.36(2) 


0.492(12) 


0.592(1) 
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